GroundwaterUpdate Subroutine

public subroutine GroundwaterUpdate(time, hillslopeFlux, percolation, capflux)

update groundwater head with Darcy law and mass conservation in two dimensions, with a macrocopic cellular automata approach

            +-----+
            |  N  |
      +-----+-----+-----+
      |  W  |  C  |  E  |
      +-----+-----+-----+
            |  S  |
            +-----+

References: Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic cellular automata for groundwater modelling: a first approach. Environmental modelling & software, 26(5), 634–643.

Arguments

Type IntentOptional Attributes Name
type(DateTime), intent(in) :: time
type(grid_real), intent(in) :: hillslopeFlux

flux from hillslope (m3/s)
type(grid_real), intent(in) :: percolation

depp infiltration from soil balance (m/s)
type(grid_real), intent(in) :: capflux

capillary flux from soil balance (m/s)

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: areaOfCell

cell area (m2)
real(kind=float), public :: fluxEC

flux East-Centre (m3/s)
real(kind=float), public :: fluxNC

flux North-Centre (m3/s)
real(kind=float), public :: fluxNESW

total flux from the 4 directions (m3/s)
real(kind=float), public :: fluxSC

flux South-Centre (m3/s)
real(kind=float), public :: fluxWC

flux West-Centre (m3/s)
real(kind=float), public :: gToR

groundwater to river discharge (m3/s)
integer(kind=short), public :: i
integer(kind=short), public :: j
integer(kind=short), public :: k
real(kind=float), public :: rToG

river to groundwater discharge (m3/s)
real(kind=float), public :: rechargeQ

vertical recharge (m3/s)
real(kind=float), public :: transEC

mean transmissivity East-Centre (m2/s)
real(kind=float), public :: transNC

mean transmissivity North-Centre (m2/s)
real(kind=float), public :: transSC

mean transmissivity South-Centre (m2/s)
real(kind=float), public :: transWC

mean transmissivity West-Centre (m2/s)

Source Code

SUBROUTINE GroundwaterUpdate   & 
  !
  ( time, hillslopeFlux, percolation, capflux )       

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime),  INTENT(IN) :: time
TYPE (grid_real), INTENT(IN) :: hillslopeFlux !!flux from hillslope (m3/s)
TYPE (grid_real), INTENT(IN) :: percolation !! depp infiltration from soil balance (m/s)
TYPE (grid_real), INTENT(IN) :: capflux !! capillary flux from soil balance (m/s)


!local declarations: 
INTEGER (KIND = short) :: i, j, k
REAL (KIND = float) :: rechargeQ !!vertical recharge (m3/s)
REAL (KIND = float) :: transNC !!mean transmissivity North-Centre (m2/s)
REAL (KIND = float) :: transEC !!mean transmissivity East-Centre (m2/s)
REAL (KIND = float) :: transSC !!mean transmissivity South-Centre (m2/s)
REAL (KIND = float) :: transWC !!mean transmissivity West-Centre (m2/s)
REAL (KIND = float) :: fluxNC  !!flux North-Centre (m3/s)
REAL (KIND = float) :: fluxEC  !!flux East-Centre (m3/s)
REAL (KIND = float) :: fluxSC  !!flux South-Centre (m3/s)
REAL (KIND = float) :: fluxWC  !!flux West-Centre (m3/s)
REAL (KIND = float) :: fluxNESW !! total flux from the 4 directions (m3/s)
REAL (KIND = float) :: areaOfCell !! cell area (m2)
REAL (KIND = float) :: rToG !!river to groundwater discharge (m3/s)
REAL (KIND = float) :: gToR !!groundwater to river discharge (m3/s)

!-------------------------------------end of declarations----------------------

!reset counter variables
volumeNeumann = 0.
volumeDirichlet = 0.
volumeRecharge = 0.
volumeChange = 0.
volumeResidual = 0.

!update boundary condition
CALL GroundwaterParameterUpdate (time)


!save head1 into head0
DO k = 1, basin % nAquifers
    basin % aquifer (k) % head0 =  basin % aquifer (k) % head1
END DO

!update flux from hillslope of first (surface) aquifer
DO i = 1, basin % aquifer (1) % domainBC % idim
    DO j = 1, basin % aquifer (1) % domainBC % jdim
        IF ( basin % aquifer (1) % domainBC % mat (i,j) == &
             BC_NEUMANN ) THEN
            
            basin % aquifer (1) % valueBC % mat (i,j) = &
              hillslopeFlux % mat (i,j)
            
        END IF
    END DO
END DO

!update aquifers head
DO k = 1, basin % nAquifers
    
    !compute local transmissivity
    DO i = 1, basin % aquifer (k) % domainBC % idim
      DO j = 1, basin % aquifer (k) % domainBC % jdim
          IF ( basin % aquifer (k) % domainBC % mat(i,j) /=  &
			     basin % aquifer (k) % domainBC % nodata ) THEN
              
             !transmissivity: for confined aquifer = ks * aquifer depth
             transmissivity % mat (i,j) = MIN ( &
                 basin % aquifer (k) % Ks % mat (i,j) * &
                ( basin % aquifer (k) % head0 % mat (i,j) - &
                  basin % aquifer (k) % bottom % mat (i,j) ) , &

                  basin % aquifer (k) % Ks % mat (i,j) * &
                ( basin % aquifer (k) % top % mat (i,j) - &
                  basin % aquifer (k) % bottom % mat (i,j) )   &
								               )
                 
                 IF ( transmissivity % mat (i,j) <= 0. ) THEN
                     transmissivity % mat (i,j) = 0.01
                 END IF
          END IF
      END DO
    END DO
    
    
    !update head1
    DO i = 1, basin % aquifer (k) % domainBC % idim
      DO j = 1, basin % aquifer (k) % domainBC % jdim
          IF ( basin % aquifer (k) % domainBC % mat(i,j) == &
               ACTIVE_CELL ) THEN
              
             !harmonic mean transmissivity
             transNC = 2.0 * transmissivity % mat(i-1,j) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity  %mat(i-1,j) + &
                             transmissivity  %mat(i,j) )
             transSC = 2.0 * transmissivity % mat(i+1,j) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i+1,j) + &
                             transmissivity % mat(i,j) )
             transEC = 2.0 * transmissivity % mat(i,j+1) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i,j+1) + &
                             transmissivity % mat(i,j) )
             transWC = 2.0 * transmissivity % mat(i,j-1) * &
                             transmissivity % mat(i,j) / &
			               ( transmissivity % mat(i,j-1) + &
                             transmissivity % mat(i,j) )
             
             !flux from North
             IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == &
                  BC_NEUMANN ) THEN
                 fluxNC = basin % aquifer (k) % valueBC % mat(i-1,j)
                 volumeNeumann = volumeNeumann + fluxNC * dtGroundwater
             ELSE
                 fluxNC = transNC * &
                     ( basin % aquifer (k) % head0 % mat (i-1,j) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxNC * dtGroundwater
                 END IF
             END IF
             
             !flux from East
             IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == &
                  BC_NEUMANN ) THEN
                 fluxEC = basin % aquifer (k) % valueBC % mat(i,j+1)
                 volumeNeumann = volumeNeumann + fluxEC * dtGroundwater
             ELSE
                 fluxEC = transEC * &
                     ( basin % aquifer (k) % head0 % mat (i,j+1) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxEC * dtGroundwater
                 END IF
             END IF
             
             !flux from South
             IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == &
                  BC_NEUMANN ) THEN
                 fluxSC = basin % aquifer (k) % valueBC % mat(i+1,j)
                 volumeNeumann = volumeNeumann + fluxSC * dtGroundwater
             ELSE
                 fluxSC = transSC * &
                     ( basin % aquifer (k) % head0 % mat (i+1,j) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxSC * dtGroundwater
                 END IF
             END IF
             
             !flux from West
             IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == &
                  BC_NEUMANN ) THEN
                 fluxWC = basin % aquifer (k) % valueBC % mat(i,j-1)
                 volumeNeumann = volumeNeumann + fluxWC * dtGroundwater
             ELSE
                 fluxWC = transWC * &
                     ( basin % aquifer (k) % head0 % mat (i,j-1) - &
                       basin % aquifer (k) % head0 % mat (i,j) )
                 
                 IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == &
                      BC_DIRICHLET ) THEN
                     volumeDirichlet = volumeDirichlet + &
                                       fluxWC * dtGroundwater
                 END IF
             END IF
             
             !total flux
             fluxNESW = fluxNC + fluxEC + fluxSC + fluxWC
             
             !area of the cell
             areaOfCell = CellArea (transmissivity, i, j) 
             
             !vertical recharge
             IF ( ALLOCATED (percolation % mat ) ) THEN
                 IF ( k == 1 .AND. percolation % mat (i,j) /= &
                                   percolation % nodata ) THEN
                     rechargeQ = ( percolation % mat (i,j) - &
                                   capflux % mat (i,j) ) * areaOfCell
                     volumeRecharge = volumeRecharge + &
                                      rechargeQ * dtGroundwater
                 ELSE
                     rechargeQ = 0.
                 END IF
             ELSE
                 rechargeQ = 0.
             END IF
             
             !upward flux from lower aquifer
             IF (k /= basin % nAquifers ) THEN
			    fluxUpward % mat(i,j) = &
                    basin % aquifer (k) % ksAquitard % mat(i,j) * &
                ( ( basin % aquifer (k+1) % head0 % mat(i,j) -  &
                    basin % aquifer (k) % head0 % mat(i,j) ) / &
                  ( basin % aquifer (k) % bottom % mat(i,j) - &
                    basin % aquifer (k+1) % top % mat(i,j)  ) )
             ELSE
			    fluxUpward % mat (i,j) = 0.0
             END IF
             
             !update head1  
             IF ( k == 1) THEN
                 fluxDownward % mat (i,j) = 0.
                 IF ( riverGroundwaterInteract ) THEN
                     rToG = riverToGroundwater % mat (i,j)
                     gToR = - groundwaterToRiver % mat (i,j)
                 ELSE
                     rToG = 0.
                     gToR = 0.
                 END IF
             ELSE
                 rToG = 0.
                 gToR = 0.
             END IF
             
             basin % aquifer (k) % head1 % mat(i,j) = &
                 basin % aquifer (k) % head0 % mat(i,j) + &
                 ( fluxNESW + rechargeQ + rToG + gToR ) * &
                 dtGroundwater / areaOfCell / &
                 basin % aquifer (k) % sy % mat(i,j) + &
                 ( fluxUpward % mat(i,j) + fluxDownward % mat(i,j) ) * &
                 dtGroundwater / basin % aquifer (k) % sy % mat(i,j)
             
             IF ( ISNAN ( basin % aquifer (k) % head1 % mat(i,j) ) )  THEN
               basin % aquifer (k) % head1 % mat(i,j) = &
                   basin % aquifer (k) % head0 % mat(i,j)
                volumeChange = 0.
             ELSE
                !volume change
                 volumeChange = volumeChange + &
                 ( basin % aquifer (k) % head1 % mat(i,j) - &
                   basin % aquifer (k) % head0 % mat(i,j) ) * &
                   basin % aquifer (k) % sy % mat(i,j) * areaOfCell 
             END IF
             
             !swap fluxUpward and fluxDownward
             fluxDownward % mat (i,j) = - fluxUpward % mat (i,j) 
             
             !update vadose zone depth
             IF ( k == 1) THEN
                 vadoseDepth % mat (i,j) = & 
                     basin % aquifer (k) % top % mat(i,j) - &
                     basin % aquifer (k) % head1 % mat(i,j)
             END IF
             
          END IF
      END DO
    END DO
    
END DO

!compute balance error
volumeResidual = volumeRecharge + volumeDirichlet + &
                 volumeNeumann  + volumeRiverToGroundwater - &
                 volumeGroundwaterToRiver - volumeChange

!write point data
IF ( time == timePointExport ) THEN
   CALL GroundwaterPointExport ( time )
END IF


RETURN
END SUBROUTINE GroundwaterUpdate